Fast chemical reaction in two-dimensional Navier-Stokes flow: 

Initial regime 
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This paper studies an infinitely fast bimolecular chemical reaction in a two-dimensional bi-periodic 
Navier-Stokes flow. The reactants in stoichiometric quantities are initially segregated by infinite 
gradients. The focus is placed on the initial stage of the reaction characterized by a well-defined 
one dimensional material contact line between the reactants. Particular attention is given to the 
effect of the diffusion k of the reactants. This study is an idealized framework for isentropic mixing 
in the lower stratosphere and is motivated by the need to better understand the effect of resolution 
on stratospheric chemistry in Climate-Chemistry Models. 

Adopting a Lagrangian straining theory approach, we relate theoretically the ensemble mean of 
the length of the contact line, of the gradients along it and of the modulus of the time derivative of 
the space-average reactant concentrations (here called the chemical speed) to the joint probability 
density function of the finite time Lyapunov exponent A with two times r and r. The time j 
measures the stretching time scale of a Lagrangian parcel on a chaotic orbit up to a finite time 
t, while r measures it in the recent past before t, and r in the early part of the trajectory. We 
show that the chemical speed scales like and that its time evolution is determined by rare large 
events in the finite time Lyapunov exponent distribution. The case of smooth initial gradients is 
also discussed. The theoretical results are tested with an ensemble of direct numerical simulations 
(DNS) using a pseudospectral model. 



I. INTRODUCTION 

The stratospheric ozone chemistry resulting from 
Climate-Chemistry Models is thought to be sensitive to 
the spatial resolution. It was shown by [l[ that the sim- 
ulated spring ozone depletion inside the polar vortex is 
very sensitive to the horizontal grid size. However 0, HJ 
pointing out some flaws in the former work, suggested 
that resolution is not crucial for ozone depletion inside 
the polar vortex during cold enough winters because chlo- 
rine, the relevant catalyst for ozone destruction, is to- 
tally activated regardless the resolution. Nevertheless, 
they suggested that at the outer edge of the vortex (the 
surf zone), where mixing is important, the filamentary 
structures exhibited by chemical fields (e.g. are not 
represented by low resolution models. The deactivation 
of polar vortex chlorine by low-latitude nitrogen oxide, 
a process controlling ozone concentrations at the outer 
edge of the mid-winter Arctic polar vortex, was stud- 
ied numerically by Q. Assuming two dimensional mix- 
ing on isentropes on time scales smaller than two weeks 
and using reanalysis data to advect chemicals, they found 
that the production of chlorine was strongly dependent 
on the tracer diffusion coefficient. They proposed that 
the product's concentration scales like k p ^ where p(t) is 
a positive decreasing function of time which depends on 
initial conditions. This problem was addressed from a 
theoretical point of view by @ which showed, for an in- 
finitely fast bimolecular chemistry, that the function p(t) 
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is given by 1 — D(t)/2 where D(t) is the box counting 
fractal dimension of the contact line between the reac- 
tants, defined as the zero isoline of a tracer <fi equal to 
the difference between the two reactants' fields. Their 
main assumption on the geometric configuration of <fi is 
that of an on/off field, which allows to link the slope of 
the tracers' variance spectrum to the box counting frac- 
tal dimension of the contact line, and the variance to the 
first moment of the modulus of </>. This interesting ap- 
proach is however limited by the lack of realism of the 
on/off fields assumption. 

Here, in the absence of this assumption, we propose 
to focus on the case where the contact line is a mate- 
rial line unaffected by diffusion (fractal dimension equals 
to one). This is true during the early stage of the re- 
action, before tracer filaments start to merge under the 
action of diffusion. To our knowledge, a detailed anal- 
ysis of this regime has not appeared in the literature 
despite its relevance to the atmosphere on time scales 
of several days to weeks. We develop a mathematical 
framework which relates the effect of diffusion on the re- 
actant concentration and its time evolution to the statis- 
tics of Lagrangian straining properties (LSP) of advected 
parcels in the flow. This approach has been widely used 
to describe the asymptotic decay of passive tracers in 
the Batchelor regime of turbulence or in chaotic advec- 
tion (for Lagrangian straining theories and further de- 
velopments, see |7H13|). In addition, this approach was 
recently applied to the long term decay of fast reacting 
chemicals by [i~H ]. 

Our assumptions are those of a two dimensional statis- 
tically isotropic, homogeneous and stationary non linear 
Navier-Stokes flow. We use ensembles of direct numer- 
ical simulations, each one corresponding to a different 
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initial condition on the vorticity, to verify the analyti- 
cal relations between the LSP and chemistry. Although 
this flow gives a very simplified representation of strato- 
spheric mixing, it can be argued that it is relevant for 
scales larger than approximately 40km ((la]). We vary 
the diffusion coefficient K of tracers to study the effect of 
resolution, employing a similar approach to Q . This 
is justified by noting that the smallest scales of the flow 
are determined by the balance of advective and diffu- 
sive processes and thus scale like . Considering that 
small-scale tracer structures are generated by the large- 
scale field, the viscosity of the field is chosen larger than 
the diffusion k. Hence, tracers evolve in a smooth veloc- 
ity field, which allows to differentiate it at the tracers' 
small scales and interpret their behavior in the frame- 
work of Lagrangian chaos. It has been shown in [HI Il7j 
that the concept of chaotic advection, where a spatially 
coarse flow produces chaotic tracer trajectories, was ap- 
plicable to two-dimensional mixing in the stratospheric 
surf zone. In addition it has been argued ([H|) that in 
barotropic, beta-plan two-dimensional turbulence, rela- 
tively coarse velocity fields reproduce quite accurately 
the fine structures of the tracer field when the spectrum 
of energy is steeper than fc -3 , which is relevant both in 
the stratosphere ([111) and in the enstrophy cascade in 
two-dimensional turbulence ( [20j| ) . 

We focus on the initial regime of an infinitely fast 
chemical reaction between two segregated rcactants in 
stoichiometric quantities. The main emphasis is placed 
on the case where the reactants are initially separated by 
a sharp gradient, while the case of a smooth gradient is 
briefly discussed. Figure Q] illustrates this regime. With 
T the integral time scale of the flow, the contact line does 
not depend on diffusion at 4 = 1 and ^ = 3, but gra- 
dients become clearly smoother when diffusion increases. 
The time span of this regime depends on the diffusion: at 
|t = 8 the contact line seems to be the same for Prandtl 
numbers Pr = 16 and Pr = 128 but is clearly different 
for Pr = 1. The Prandtl number is defined as the ratio of 
the tracer diffusion k with the fluid viscosity v and when 
the diffusion is larger, filaments merge earlier, making the 
contact line dependent on diffusion at a smaller time. 

This approach is relevant to the chlorine deactivation 
at the outer edge of the winter time polar vortex, which is 
very fast compared to advective and diffusive time scales 
(Q). It is also of general interest in isolating and investi- 
gating the effect of two-dimensional turbulent mixing on 
chemical reactions. A separate paper, in preparation, will 
focus on the case of a more complicated contact line (box 
counting fractal dimension between 1 and 2), which corre- 
sponds to the intermediate and time asymptotic regime. 

This paper is organized as follows. Section II. de- 
scribes our approach and methodology. We show that 
with infinitely fast chemistry, average concentrations of 
reactants and product are simple linear functions of the 
first moment of the modulus of the passive tracer concen- 
tration <p defined as the difference between the reactant 
fields. This approach is rather general in the study of 



infinitely fast bimolccular reactions ( [y, [Tj, [2l|, [22[ ) . In 
particular, this implies that the reaction is controlled by 
the diffusive flux across the isoline = 0, noted C. The 
importance of the behavior of C for chemistry in com- 
plex flows or complex geometric configurations of chem- 
ical fields has been highlighted in [2^, [H| respectively. 
This section also describes the numerical model and the 
simulated flow, including the spatial configuration and 
the probability density function (pdf) of the FTLE. Sec- 
tion III. describes the theoretical and numerical results. 
We derive analytical expressions for the lengthening of 
£, for the gradients advected along C and finally for the 
diffusive flux across L, the latter being equal to the time 
derivative of the space-average reactants concentrations. 
We compare the theory to ensembles of numerical sim- 
ulations. Finally, conclusions are drawn in section IV. 



II. METHODOLOGY 

A. The limit of infinite chemistry 

We consider the bimolccular chemical reaction A + 
B — > C with C A , Cb and Cc, the concentrations of A, 
B and C respectively. Eulerian equations describing the 



evolution of Ci(x, t), 
are: 



A, B, C, in the flow u = («, v) 



dC A 
dt 
dC B 
dt 
dC c 
dt 



u • VC A = kV 2 C a - k c C A C B 



- u • VC B = kV 2 C b - k c C A C B 
u • VC C = ^ 2 C C + k c C A C B , 



(fa) 
(lb) 
(lc) 



where k c is the chemical reaction rate and k the diffu- 
sion, which is assumed equal for all tracers. The quantity 
<p = C' A — C B is a passive tracer which obeys the simple 
advection-diffusion equation 



^ + u • V</» = kV 2 <; 
at 



(2) 



We assume that </> has zero spatial average, which is 
equivalent to having the rcactants in stoichiometric bal- 
anced ratio. Under the fast chemistry hypothesis (k c — > 
oo), we can assume without loss of generality that the 
reactants A and B are segregated (i.e. A and B do not 
overlap spatially). In fact, even if they are co- located at 
time t = 0, they can not coexist at a later time t > since 
they react instantaneously where both fields are together 
non-zero. It follows that: 



C A (x,t) = ct)(x,t) andC B (x,i)=0 if <f>(x, t) > 
C B (x,t) = -<j){x,t) andC A (x,i)=0 if<£(x,i)<0 

(3) 
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FIG. 1. (Color online) Reactant fields in a bi-periodic domain [— tv,tt] 2 . Red (or light grey), positive values, and blue (or dark 
grey), negative values, show the two reactants A and B. From left to right Pr = 1, 16, 128 and from top to bottom ^ = 1, 3, 8. 
The Prandtl number Pr is defined as the ratio of the viscosity of the fluid to the tracer diffusion. Since the viscosity is fixed, 
an increasing Pr means a decreasing diffusion. T is the integral time scale of the flow. 



Defining with an over-bar the average over the whole do- 
main, we have: 



C a = C b = y 



Cr 



m=$)\-\<p\ 



(4a) 
(4b) 



Consequently, studying the decay of the reactants of 
an infinitely fast chemical reaction in stoichiometric bal- 
anced ratio is equivalent to studying the decay of the first 
moment of the modulus of a passive tracer <f> of zero spa- 
tial average. For an incompressible flow, it can be shown 
with the divergence theorem that the decay rate of the 
total reactant quantity for an infinite reaction equals the 
diffusive flux across the contact line between A and B, 
namely C = {x|0(x) = 0}, oriented in a counterclockwise 



direction around the area where A is located: 

dte l .dW\ 



A 



dC A 
dt 



A 



dt 



2 dt 



ndl 



C(t) 



cm 



|V0|dZ, (5) 



where A is the total area of the domain and n the vector 
normal to C and pointing outside the area where A is 
located. The contact line is by definition an isoline of 4>, 

which gives the last equality in ([5]). Hereafter, — is 
called the chemical speed. 



B. The numerical model 



1. The flow 



The numerical model integrates the vorticity equation: 



~dt 



u • Vw = F - R uj + vV 2 w 



(6) 
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where ui = V x u is the vorticity, F the forcing term, Rq 
the Rayleigh friction and v the viscosity. The equation is 
integrated in a bi-periodic domain (x, y) € [— 7r,7r] on a 
512x512 grid using the pseudo-spectral method. The fast 
Fourier transforms are provided by FFTW (UK). The 
Fourier series are truncated at K max = 512/3 to avoid 
aliasing. The time stepping algorithms are leap-frog for 
the advection and Crank- Nicholson for the viscosity. The 
computational mode is dissipated by a weak Robert fil- 
ter with parameter 0.001. The forcing term F has the 
following form in Fourier space: 



Fu = 



0.002 if k = (±3, 0) and k = (0, ±3) 
otherwise 



(7) 



The energy tends to concentrate in the largest scales of 
the flow because of the inverse energy cascade inherent to 
two-dimensional turbulence. As a consequence, we use a 
Rayleigh friction term with Rq = 0.0002 in the vorticity 
equation ([B]) to balance the injection of energy through F. 
The viscosity is v ~ 5.57 x 10 -4 and results in a Reynolds 
number Re of the order of 10 3 . It has deliberately been 
chosen to be relatively low for reasons explained in the 
Introduction. 

A snapshot of the vorticity field is depicted in fig- 
ure H] (top left). With brackets for an ensemble av- 
erage, the flow has an RMS velocity (u • u) 2 ~ 0.08 



and a mean enstrophy Z 



0.009, which cor- 



responds to an advective time scale T = Z~? ~ 10 

(|18|). Hereafter, T is used to normalize time and 

can also be estimated from the mean strain rate (S), 

1 

where S = ) 2 + \ + f^'V (here expressed 

in Cartesian coordinates for an incompressible flow). In 
two dimensional turbulence, we have ((2S) 2 ) = (lu 2 ) (e.g. 
[2(| ) . The distribution of the strain is close to a Rayleigh 
distribution Q, as a consequence we have (S 2 ) w ^(S) 2 . 
Finally, we have T w x/JjtW' ^he term \/l being of 
the order of unity, the turnover time can be evaluated 
from ^7=r. The mean strain rate is about 0.05 in our 

flow (see figure [2] for the whole distribution), which gives 
approximately the same estimate as Z~? for T. 



2. Finite Time Lyapunov Exponents (FTLE) 

a. Definition and properties The finite time Lya- 
punov exponent is defined as the rate of exponential in- 
crease of the distance between the trajectories of two fluid 
parcels that are initially infinitely close. If Sl(t) is the dis- 
tance at time t between two parcels that start at x and 



1 It would be exactly a Rayleigh distribution if the velocity deriva- 
tives had Gaussian statistics and where statistically independent. 



x + SIq at time t — 0, then the FTLE X(x,t) at x over 
the time interval t is 

\(x, t) = - max lim { In if^U , (8) 

t a |M |-)V I (H 1 

where the maximum is calculated over all the possible 
orientations a of SIq. The unit vector with the orienta- 
tion ^ + (x, i) of SIq at the maximum is called a singular 
vector and we note it i/)+(x.,t) = (cosip+,smijj+). It de- 
fines a Lagrangian straining direction. It follows from 
© that the FTLE converges to the strain rate as t — > 0. 
For large times, the large deviation theory suggests that 
the FTLE pdf P\ in chaotic flows without KAM (Kol- 
mogorov, Arnold, and Moser) surfaces ((27[) can be well 
approximated by: 



P\(t, A) 



*G"(A ) 



2n 



expHG(A)), 



(9) 



where G(A), the Cramer or rate function, is concave 
with its minimum at Ao satisfying G(Ao) = G'(Ao) = 0. 
Moreover, Ao is the infinite-time Lyapunov exponent: 
lim P\(t, A) = 6(Xq — A) where 5 is the Dirac delta func- 

t— >oo 

tion. The convergence of the Lypunov exponents is very 
slow and typically algebraic in time ([H|). The form © 
has been numerically verified and is widely used to ap- 
proximate the asymptotic form of FTLE pdfs in simple 
ergodic flows with chaotic advection (e.g 0, EH EI ) ■ 

b. Computation and description The distance dl be- 
tween two trajectories initially infinitely close is solution 
of 



^ - S.0I = 0, 

dt 



(10) 



where the tensor S = Vu(X, t) is the velocity gradient 
tensor along a trajectory X(x, t). The distance dl can 
be calculated by 61 — M.SIq, where the resolvent ma- 
trix M is solution of — SM = and is equal to the 
identity at t = 0. The finite time Lyapunov exponent A 
is times the log of the largest eigenvalue of [M T M], 
with the singular vector ip + the associated eigenvector. 
The FTLE are obtained using the method described in 
[29l ] from the trajectories computed offline using a fourth 
order Rungc-Kutta scheme with a trilinear interpolation 
on the velocity field. The time step is 0.1, which corre- 
sponds to a hundredth of the turnover time. The tensor 
S is calculated along the trajectories to obtain M and 
consequently A and 

We estimate the FTLE pdfs as normalized histograms 
over 100 realizations of the flow, differing by their ini- 
tial vorticity field. We initialise a trajectory at every 
grid point of our 512 2 domain, which results of a total of 
about 26 x 10 6 trajectories calculated. Each realization 
is run for a time span of 25T. The FTLE pdfs are shown 
at different times in figure [2] The variance of the FTLE 
decreases with time while the peak of the distribution 
converges toward X max ~ 0.02. The FTLE are signifi- 
cantly smaller than the strain rates probably because of 
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FIG. 2. (Color online) Density P\ of the finite-time Lyapunov 
exponents shown at different times between t — and t = 
25T. As the time increases, the density shifts toward smaller 
values. 



vorticity that inhibits the stretching of material elements 
and because of the reorientation of the local strain 
axis along a trajectory (e.g. |3lj|). In order to estimate 
whether these pdfs are asymptotically well approximated 
by ©, we define 



Ge(X,t) 



ln(P A (t,A)) 



t 



hit 
~2t 



A e (t) 



t 



(11) 



where A e (t) is chosen such that minG e (A,t) = 0. Fig- 
ure [3] shows the time evolution of G e . The convergence 
for large values of A, typically larger than the mean (A) 
is satisfactory. However, the convergence for small val- 
ues is much slower. It is particularly difficult to get the 
Cramer function for small values of A ([IH). This is not 
a concern for the present study because only values of A 
larger than their average are relevant. Nevertheless, we 
can get an estimate of the Cramer function assuming it 
is symmetric, as represented on figure [3] We have fitted 
the average of the Cramer function for times larger that 
15T for values of A larger than 0.02 with a second order 
polynomial (gaussian approximation) and obtained the 
estimate for values smaller than (A) using symmetry. 

The FTLE maps are shown in figure|4] For small times 
the strain field is dominated by large scales because two 
dimensional dynamics do not allow a cascade of energy 
toward smaller scales. However filamentary structures 
appear shortly, becoming finer and finer until they reach 
the resolution of the Eulerian model (the trajectories are 
initiated at every grid point of the Eulerian model). It 
is interesting to note that we get very similar structures 
as [26}, despite our much coarser velocity field. In par- 
ticular, in the direction perpendicular to the FTLE fila- 



FIG. 3. (Color online) Function G e (A,t) plotted at different 
times (0 < t/T < 25). G e is defined such that Pa - plotted in 
figure[2]- can be written oc — tG e (A, i) with minG e = 0. As the 

time increases, the function G e shifts toward smaller values. 
We note the asymmetry of G e and the faster convergence for 
FTLE larger than their ensemble mean. The time asymptotic 
form of G e is the Cramer function G corresponding to the 
longtime FTLE pdf P\. An estimate of G is given by the 
blue circles, using a method detailed in the text. 



ments, it can be observed than initially close trajectories 
can have totally different FTLE. This might be a mani- 
festation of chaotic advection. 

It has been argued, in ergodic systems, that the singu- 
lar vectors converge exponentially in time ([Hj]), faster 
than the Lyapunov exponents, whose convergence is al- 
gebraic ([1^, |28|, HH). The "freezing" of the large scale 
patterns in the FTLE maps (figure @| may be interpreted 
as a manifestation of the convergence of the singular vec- 
tors. In fact, [28j argued that in ergodic and conser- 
vative chaotic dynamical systems, the Lyapunov expo- 
nents varies slowly along lines (the s lines) which defines 
the stable direction in which neighboring points asymp- 
totically converge. The filamentary structures in figure 
H] may be interpreted as being these s lines. This has 
been verified experimentally through the computation of 
the singular vectors (not shown), their convergence be- 
ing particularly fast for trajectories originating in areas 
of the flow dominated by strain. 

In the theoretical developments of part III.B and III.C, 
we will neglect the time evolution of the Lyapunov vec- 
tors and will only take into account the time evolution of 
the Lyapunov exponents. 




FIG. 4. (Color online) Maps of FTLE calculated at different times and displayed at starting locations of trajectories in the 
bi-periodic domain [— tv,tv] 2 . In the top row are shown the vorticity (left) and the strain (right) at t = 0. In the following 
panels, ordered from left to right and top to bottom are shown the FTLE maps for ^ = 1, 2, 4, 8, 16, 24. 
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3. The tracers 

The passive scalar <f> is integrated with equation us- 
ing the same numerical scheme as for the vorticity. The 
numerical simulations are performed for eight different 
Prandtl numbers Pr = * = 2 l for < i < 7. Con- 
sequently the Peclet number Pe = PrRe, which mea- 
sures the ratio of the advective to the diffusive time scale, 
ranges from 10 3 to 10 5 . 

We use two different initial conditions on the tracer for 
(x,y) G [-7T, ?r] 2 : 

4>{x,y,t = 0) = A sgax 

for infinite initial gradients (12) 



(x, y, t = 0) = —r cos x cos y 

for smooth initial gradients 



(13) 



where sgnx is the sign of x and A$ is twice the initial 
domain average concentration of both A and B in the 
box. The first initial condition allows represents the case 
of sharp (actually infinite) gradients separating areas 
of well mixed reactants and the second one the case of 
smooth gradients. 



C. The ensemble analysis 

For each value of Prandtl number and for each ini- 
tial condition, we run an ensemble of 34 simulations (or 
members). Each member is defined by different initial 
condition on the vorticity, taken from a long simulation 
of the statistically stationary flow solution of ©. 



III. THEORETICAL AND NUMERICAL 
RESULTS 

Our goal here is to describe and understand the initial 
evolution of the first moment of \<fi\. In other words, we 
would like to integrate ([5]). We first consider how a ma- 
terial line stretches in a Lagrangian framework (III. A), 
and then how gradients on the contact line evolve under 
the action of both the diffusion and the flow along a La- 
grangian trajectory (III.B). Paragraph III.C deals with 
the chemical speed. We focus on the initial condition 
where the reactants are separated by a sharp gradient 
before discussing the case of smoother gradients (III.D). 



A. Lengthening of the contact line C 

1. Theory 

We consider a line element <51o along the contact line 
C(t = 0) = Co- Its coordinates are (%o(cosa,sina). The 



angle a is the initial orientation of the line element. It is 
transformed at time t into an element SI = M.61q whose 
norm is: 

\Sl\ = [8lTM T M5lo] ' 



= \Sl \ [e 2At cos 2 (^+ - a) + e~ 2Xt sin 2 (?/>+ - a)} 5 . 

(14) 

The resolvent matrix M was introduced in II.B.2.b. The 
angle a and consequently the angle 7 = — a between 
the initial orientation and the singular vector can be as- 
sumed uniformly distributed between and 2n and sta- 
tistically independent of the chaotic orbit because the 
contact line is chosen arbitrarily with respect to the flow. 
Integrating over the Lyapunov exponent A, the angle 7 
and the initial contact line, gives the ensemble average 
(L) of the length L of C (brackets are for ensemble aver- 
ages). With P\ the probability density distribution of A, 
we have: 



POO /*2lT 

(L) =L / [, 

J\=0 ^7=0 



e 2At cos 2 7 + e- 2At sin 2 7 ] 5 



P x (t 7 X)dX^. (15) 



The length Lq is the initial length of the contact line. 
Equation (|15|) gives the actual length with no diffusion. 
Given the chaotic and closed (periodic) nature of the 
flow, we can only neglect diffusion as long as the contact 
line has not folded on itself. Indeed, when two filaments 
of C are brought together at a distance smaller than the 
diffusive cutoff, they merge under the action of diffusion. 
The time span of the regime where (|15[) is expected to be 
valid can be approximated with the mix-down time T m i X 
from the the largest scale L of the flow to the diffusive 
cutoff L K which is, according to 34|, \ \a(L/L K ), where 
A is the thinning rate of a fluid clement, i.e the Lyapunov 
exponent. It follows that T m i X depends on the trajectory 
we are considering. To obtain an estimate of T m i X , we 
use A sa (S) and L K ? 



TlnPe = TlnRePr, 



(16) 



The length (L) can be approximated by when 
we neglect the sine term in (|15[) . i.e. when the con- 
tact line elements have equilibrated with the flow: their 
length converge to a function that grows exponentially 
at a rate given by the FTLE, the initial orientation a 
of the contact line being "forgotten" . This is valid for 
t > 



i_ 

4<S> 



T 
2 ' 



poo p2n 

(L) ~ Le = Lq / P A (;,A)|cos 7 |e At dA 
*>t Jx=oJ 7 =o n 



At « dl 



2i n 



P x (t,X)e M d\. 



(17) 



If wc assume P\(t,X) cx e G <:( A ^)* ; with G e a concave 
positive function, integrating (flT)) with the steepest de- 



8 



scent method, we obtain: 

/>oo 

Leo, e [A-G e (A, t )] tdA = e ma XA [A- Gc (A,t)]^ (lg) 

Jo 

Asymptotically, we have: 

L E ^e Xlt , (19) 

where 

Ai = max[A - G(A)] (20) 

is the Legendrc transform of G evaluated with argument 
one. The value of Ai from our numerical estimate of G 
(figure [3]) is 0.027. 

2. Numerical results 

The theoretical predictions (L) and Le are compared 
to the numerical calculations in figure [5] The integration 
of (|15p using our numerical estimate of P\ reproduces 
very accurately the initial lengthening of the contact line 
for t < T m i X (Pr). The derivative of the mean length of 
the contact line at t = is because contracting line 
elements statistically compensate with stretching line el- 
ements due to randomness of a. The inflection of ln(L) 
around ^ = 2.5 is due to two opposite effects: the equi- 
libration of the contact line with the flow accelerates the 
growth of the line, while the shift of the FTLE pdf to- 
ward smaller values decelerates it, as shown by the Le 
curve. 

As seen on figure [5l (L) and L e have a behavior very 
close to an exponential increase at the rate Ai « 0.027 
after a couple of turnover times. This is consistent with 
the fast convergence of G e for large FTLE (figure [3]). 
A behavior close to this exponential increase can actu- 
ally be seen in the simulations with large Prandtl num- 
bers for a window of turnover times from around 4T to 
6T. Note that numerical simulations with even larger 
Prandtl numbers would have increased this time window 
only marginally since dividing the diffusion by two ex- 
tends its time span by only half a turnover time (|16p . 
The reason is that the convergence of G e , at least for 
larger than average FTLE, has a time scale close to the 
advective time scale. Further investigations are needed 
to explain this fact. 

It is worth noting that the lengthening of a material 
contour is determined by rare events in the tail of the 
FTLE distribution. The maximum max[A — G e (X,t = 

NT)] is achieved by values of A in the 42% quantile of 
the distribution for N = 2, 27% for N = 4, 13% for 
N = 7, and 3% for N = 15. Even though those events 
become exponentially rare because of the convergence of 
the FTLE pdfs toward a Dirac distribution, their con- 
tributions to the ensemble average of the contact line 
become exponentially important in the average of expo- 
nentials (|15p. 




2 4 6 8 10 

Time in units of turnover time T 

FIG. 5. (Color online) Ensemble average of the length of 
the contact line (infinite initial gradient case). The symbols 
correspond to ensemble averages of DNS for different Prandtl 
numbers Pr = 1,2,4,8,16,32,64 and 128. The black solid 
line corresponds to (L) , as estimated from the FTLE pdf using 
(|15[) . and the bold blue dashed line to Le as estimated using 
the FTLE pdf with simplified expression (jTTJ) . The light red 
dot-dashed line corresponds to an exponential increase at a 
rate Ai = max[A — G(A)] ~ 0.027, which corresponds to the 

asymptotic behavior in the inviscid limit (|19[) . The latter has 
been shifted vertically for clarity. We note the log scale in the 
y-axis. 



B. Lagrangian advection of the gradients along the 
contact line C 

In this section, we calculate the time evolution of the 
gradient of (f> along a Lagrangian trajectory on the con- 
tact line for infinite initial gradients. We take into ac- 
count the time evolution of the Lyapunov exponents. The 
singular vectors are taken equal to the forward Lyapunov 
vectors: ip+(x,t) = \&+(x). As noted in part II. B. 2, since 
the singular vectors converge rapidly as t increases, we 
expect this approximation to yield an accurate estimate 
of the gradients because the singular vectors converge 
rapidly, as noted in part II. B. 2. 

1. Advection- diffusion equation in a co-moving frame 

We consider a fluid clement on the contact line and we 
note X£ its trajectory: 

^£ = u(X £) i) with X £ (t = 0) = X 0£ , (21) 

where X £ is the initial location of the contact line el- 
ement we are following. We define a new coordinate r 
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corresponding to a frame co- moving with X £ : 



using the rescalcd time 



(22) 



e 



\re 2M cos 2 7 + t sin 2 7] 



(29) 



Writingthe concentration field X (r, t) = 0(x, t), we can We have reintroduced 7 = - a, a random and uni- 
show ( 35]) using © and dH]) that: 



^ + [u(X £ + r, t) - u(X £ , *)] • V X = ^V 2 X . (23) 

Assuming a separation of scale between velocity and 
tracer scales, we can write [u(X £ + r,t) — u(X £ ,t)] at 
the first order in r. We basically assume that the char- 
acteristic width of the contact zone is much smaller than 
the velocity scale. 



dt 



v 1 .Vu(X £ ,t)-Vx = ^V 2 X . 



(24) 



Locally, along the contact line, the concentration of <f> 
only varies in the direction perpendicular to the contact 
line, assuming that, for t < T ril i X , the curvature of the 
contact line is much larger than the width of the contact 
zone where the gradients are concentrated. This is rele- 
vant because the stirring in chaotic advection produces 
elongated structures by nature. As a consequence, as 
noted previously in a similar case ([111), the field x has 
to be of the form: 



= X(k- r,t) = x(r),t), 



(25) 



where k is a vector perpendicular to the contact line and 
■q a coordinate along k. Substituting (|2"5j) into (|2~4"]) (with 
S(i) = Vu(X £ ,i)) and equaling the zero and first order 
terms in r, we can show that f [361]): 



dk 
~dt 



S J .k = 

d 2 x 



^- K \k\ 2 - 
dt~ K |k| dlf 



(26a) 
(26b) 



Equation (|26ap is actually the equation of a wavenumber 
k advectcd with the trajectory X £ . Noting its similarity 
with ((TO)) , it is clear that the FTLE is also the max- 
imum exponential growth rate of a wavenumber k (or 
equivalently of a passive tracer gradient in the absence 
of diffusion). This is an alternate and classical defini- 
tion of FTLE ([26]). Considering the resolvent matrix 
N such that k = Nk , where k = k (— sin a, cos a) is 
the initial value of k, the finite time Lyapunov exponent 
A is the log of the largest eigenvalue of [N T N]2t with 
(— s'mtp+, cos-0+) the associated eigenvector. As a con- 
sequence, we have 

|fc| 2 =fcjN T Nfc 

= |fc | 2 [e 2At cos 2 (V>+ - a) + e- 2Xt sin 2 (V>+ - a)] . 

(27) 

With the assumption ip + (x,t) = \& + (x), equation (|26b|) 
can be written 



dx _ U 2 d2 X 
de - dr, 2 



(28) 



formly distributed angle between and 2n (sec III.A.l). 
The quantities r and r are two equivalent times defined 
as follows: 



e 2u\(u) du 



2tX(t) 



and t 



2uXM du. 



(30) 



The time r, introduced by [3] and called "equivalent 
time" by [37| evaluates the stretching time scale of a La- 
grangian parcel in the recent past because chaotic trajec- 
tories are characterized by positive Lyapunov exponents. 
Similarly, the equivalent time r measures the stretch- 
ing rate in the early part of the trajectory. As a conse- 
quence, we expect r and r to have the same statistics, 
to be asymptotically equivalent as t — > and to become 
independent at larger times. It has been argued ([37^) 
that the pdf of r converges to a time asymptotic form, 
which is suggested for our flow in figure [6] where we have 
plotted the pdf of - calculated together with the Lya- 
punov exponent on each Lagrangian trajectory (II. 2. B). 
The statistics of r (not shown), calculated the same way, 
are not distinguishable from these of r. 



2. Solution (infinite initial gradient case) 

The initial gradient along the contact line is infinite, 
while the reactants are well mixed in their respective 
domain with a concentration equal to Aq. As a con- 
sequence, we take: 

X(v,t = 0) = A sgnr]. (31) 

The solution of (l28l) with the initial condition (I3P1) is: 



e~ l ~ dl = A Q Erf 1 



V 



2fc n VKO' 



(32) 

The function Erf is the Gauss error function. It follows 
from fl25J| and (32): 



x( r i t) = Aq Erf 1 



2y^ 



k| /fct,, 



e 



■4„Erf(||„.r) 



with 



(33) 



(34) 



the unit vector normal to the contact line and 
G c - 



|k| /ko 







e 2At CO g^ j _j_ g- 



-2At 



sin 7 



Tg 2At CO g2 j _|_ r gm 2 ^ 



(35) 
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The norm |V</>£| of the gradient of the field on the 
contact line (where X = 0), i.e. at the location of the 
trajectory X£ characterized by the Lagrangian straining 
properties (A, r, r , 7) is 

|V^| = |V rX -n| r=0 = A 



:G C (t,\T,T,i). (36) 



3. Ensemble average of the gradient along the contact line 

To perform the ensemble average (\V(j)c\) of the mod- 
ulus of the gradient of cj> along the contact line, we intro- 
duce the joint pdf P of (A, r, r). As noted previously, the 
orientation 7 is assumed uniformly distributed between 
and 27r and independent from the random vector (A, r, r). 
If we consider a trajectory X£ of a contact line element 
61, the gradient on it is equal to |V$c| on a length \Sl\, 
defined in ([TO . As a consequence, with (|36|) and (|TO . 
we obtain: 



(|V0£|) = 



(l«l) 

(L) 



e 2At cos 2 7 



g -2At sm 2^, 



V re 2At cos 2 7 + r sin 2 7 



P(£, A, r, r) dA dr dr 



d/) 

2tt' 
(37) 



The integration is performed between and 00 for A, 
t and t and between and 2ir for 7. Hereafter, these 
bounds will be omitted. For times sufficiently large (t ^> 
^iy rj T), we neglect the sin 2 terms under the integral 

in (|37p and in the expression for (L) and we obtain, in 
the limit of a contact line equilibrated with the flow: 



(\V<f>c\) 



^=P X ,r{t,\,T)d\dT, (38) 



where P\. T is the time dependent joint pdf of A and r. 
The joint density of (A, -) is pictured on figure[7] The fre- 
quencies A and ^ are clearly dependent, especially when 
they are small, even at times much larger than the ad- 
vectivc time scale (e.g. t = 20T and t = 25T). The 
computation of the Spearman Rho correlation coefficient 
cle arly confirms this dependence Previous studies (c.g 
0, [37| ) have assumed the independence between A and 
t at times much larger than the Lagrangian correlation 
time, here shorter than or of the order of the advective 
time scale. This may be apporpriate in simple ergodic 
chaotic flows. However, two dimensional Navier-Stokcs 
flows, including two-dimensional turbulence, exhibit co- 
herent structures (vortices, filaments of vorticity, etc..) 
that seem to prevent this independence to be achieved. 
Nevertheless, the dependence is weaker for large values 
of A, which precisely dominate the integral (|38[) . Approx- 
imating P\ T by the product of its marginal densities Pa 
and P T , we obtain that |V0£| can be approximated by 
the simple expression —4= 
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FIG. 6. (Color online) Probability density function of - plot- 
ted at different times. The equivalent time r is defined in 
(|30[) . As the time increases, the density shifts toward smaller 
values. 



4- Comparison with the numerical results 

The ensemble average of the modulus of the gradient 
along the contact line have been calculated on the 34 
ensemble members and for the whole range of Prandtl 
numbers Pr = - = 2* for < i < 7. We calculate A, 
t and r on each trajectory, which permits the numerical 
integration of (|37|) and (|38|) . Numerical results are dis- 
played on figure [5] and are compared with the theoretical 
results of the previous paragraph. The joint statistics of 
(A, r, r) are referred to as the Lagrangian straining prop- 
erties (LSP). 

The ensemble averages of the gradient calculated from 
the DNS and multiplied by v™L arc shown for Prandtl 
numbers ranging from 2 to 128 in figure [8] For large 
enough diffusion (small enough Pr), the curves become 
virtually identical, showing the dependence in K~i of the 
gradient suggested by equation (|37[) . This regime seems 
to be valid for times up to 3.5T at Pr = 2 and up to 6T at 
Pr = 16. This timescale corresponds to that estimated 
by equation (fl6|) modulo a factor 2 and coincides with 
the regime where the advection alone accounts for the 
lengthening of the contact line (figure [5]). The departure 
at small diffusion comes from the fact that the infinite 
gradient hypothesis becomes inaccurate in the numerical 
simulations given the finite size of the grid. We can re- 
produce the curves at large Prandtl number by solving 
the derivative with respect to rj in (|28[) using the initial 



condition on the gradient 



dxt I 
dr] lt=o 



2.50^ 



with 8q 



a length corresponding to a grid point. We find that the 
previous developments stand with Gc (|35[) replaced by 
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IS 

FIG. 7. (Color online) Correlation between A and - as a function of time (top left) and joint pdf of (A, —), as estimated from 
the numerical simulations and plotted at different times 4 = 0.25, 2, 4, 7, 12, 20 and 25. 
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FIG. 8. (Color online) Ensemble average of the gradients 
advected with the contact line, multiplied by -^p^, in the 
sharp gradient case. The symbols correspond to ensemble av- 
erages over the 34 DNS members for different Prandtl num- 
bers Pr = 2, 4, 16, 32, 64, 128. The lines correspond to the 
calculation from the Lagrangian straining properties (LSP): 
in solid from (|37|) and in dashed from (|38|) , considering a per- 
fectly equilibrated contact line with the flow. The dot-dashed 
line is — ^ and corresponds to (|38[) with A and r statistically 
independent. We note the log scale in the time axis. 



C. Time evolution of {^r) 



1. Theory 

Having formulated the evolution of the contact line 
and the gradient, we can now express the chemical speed 

— (^jr-) by ensemble averaging ([5]). We use the expression 
of \6l\ in (fT4"]) . with ip + = for dl and the expression 
|V</> £ | in d36j for |V^|. 



v dt ' 



LqAq 



e 2Xt cos 2 7 + e~ 2Xt sin 2 7 \ 
V Te 2Xt cos 2 7 + r sin 2 7 / 

e 2Xt cos 2 7 + e~ 2Xt sin 2 7 
\J re 2Xt cos 2 7 + r sin 2 7 



P(t, A, r, r) cZA dr dr — 



(39) 



2 -^v^ // ^Px jT (t,A,r)dAdr. (40) 



The chemical speed scales like fc2 , which is a direct con- 
sequence of the scaling of the gradients like kT^ , the 
contact line length being independent of the diffusion in 
the regime considered. Indeed, comparing equation (|39[) 
with ([57)1 leads to the simple relationship between the 
ensemble means: 



-<«> = ><|v W 



(41) 



This relationship was actually previously justified and 
used to calculate the gradients (|3T|) . 



G £ , K = / / Atcos27+e ' 2At _ sin27 which is a function of 

y -^- + [re 2At cos 2 7+r sin 2 7] 

k. The expression Gc is a good approximation of Gc,k 
when the initial gradients imposed by the grid ^ are 
large compared to -^S= (t/kt can be interpreted as the 

diffusive cutoff). This is not the case for Pr = 64 and 
Pr = 128 in our simulations. 

The evolution of ^^-(\V(f>c |) estimated using (|37| with 
the LSP is provided on figure [5] It captures very well the 
behavior of ensemble mean gradients for small Prandtl 
numbers. A very slight underestimation is seen that 
could be due to numerical artifacts or to our approxi- 
mation taking the singular vectors constant in the theo- 
retical developments. We also show ^^dV^I) approx- 
imated by (|38]l . It overestimates the gradients at small 
times, a discrepancy which decreases with time as the 
contact line equilibrates with the flow. The quantity , 
also shown in figure[8]neither performs well at small times 
for the same reason as (f3"8"|) . nor at larger times because 
of the missing dependence of r with A. 



2. Numerical results 

Figure [9] shows — -^(^^) for various Prandtl numbers 
estimated from the ensemble DNS and the result of equa- 
tion (f3"9")) using LSP. Like for the gradients, the curves 
converge together when the diffusion gets larger, for times 
shorter than T m i X . The limit curve best fulfills the infi- 
nite gradient hypothesis and consequently matches very 
well the estimate from equation (|39|) calculated from the 
Lagrangian straining properties. 

The general behavior of the chemical speed can be 
interpreted in light of equation (1411) . The initial de- 
crease is mainly due to the decrease of the gradients, 
as observed previously. Then, it is dominated by the in- 
crease of the contact line, the gradients decreasing very 
slowly. Figure |H] shows the exponential increase at a rate 
Ai w 0.027. The timescale corresponding to the mini- 
mum of the chemical speed can be estimated from the 
timescale of the decrease of the gradient, which is of the 
order of tQ 



Assuming that it is the time scale for the decrease of — ;= , this es- 
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FIG. 9. (Color online) Ensemble average of the chemi- 
cal speed in the sharp gradient case divided by the dif- 
fusion y/n. The symbols correspond to numerical results 
from the 34 members ensemble, for different Prandtl num- 
bers Pr = 1,2,4,8,64,128. The solid line (calculation from 
the LSP) corresponds to (|39l) . The exponential increase at 
a rate Ai = 0.027 (dashed line) corresponds to the expected 
asymptotic regime of (J39}. 



D. Alternative initial condition on the tracers: 
smooth gradients 

The following calculations are extending the analytical 
results for an initial condition on the tracers with smooth 
gradients and are validated numerically with the initial 
condition (|13|) . We neglect the diffusion to determine the 
evolution of the gradients. In the inviscid limit, a gra- 
dient along a Lagrangian trajectory obeys the wavenum- 
ber equation (|26a[) (e.g. [2o|). whose solution is given by 
([2~T]) . Together with (|14[) . we obtain the ensemble average 
of©: 

_(M) =/3k (( Px (t,X) [e 2At cos 2 7 + e- 2At sin 2 7 ] ^dX 



(42) 



~ P K f P x (t,X)e 2Xt d\o,e m ^ 2X - G ^ t 



(43) 



,A 2 t 



timate is obtained by direct calculation of —j= from (130 1 > taking 
A rj S. The latter approximation is justified because the Lyapunov 
exponent is very close to the strain rate where the trajectory orig- 
inates for times smaller than T, as shown on figures Handll] The 
fact that we find a time scale of the order of T validates the ap- 
proximation J^Sa posteriori. 



max [2 A 

A 



(3(A)] 



(44) 



the Legendre transform of G evaluated in two and /3 = 

£o<|V0£|)(t=O) 
A 

The dependence of the chemistry on the diffusion is, 
like in the sharp gradient case, algebraic but the expo- 
nent is now 1. Our numerical simulations are consistent 
with this prediction: figure [10] shows the chemical speed 
divided by the diffusion. For small times, all the curves 
are virtually identical, which confirms the k dependence 
of the chemical speed. 

The calculation of -(^) using @2) with the pdf P\ 
reproduces very well the initial increase of the chemical 
speed, as shown on figure 1101 Interestingly, the chemi- 
cal speed has now a similar evolution as the contact line 
(figure [5]). Indeed, equations (fT5)) and (|4"2"j) are very simi- 
lar. The quantity integrated over the density of A is just 
squared in (|42j) compared to (|15|) . Using our numerical 
estimate of the Cramer function through ([44]) , A2 can be 
estimated at 0.09, which is about three times Ai. The 
chemical speed increases much faster than twice the con- 
tact line, which would be the case for a uniform Lyapunov 
exponent. In the sharp gradient case, the chemical speed 
rather scales like e Alt , because of the action of diffusion 
on the gradient. This suggests a smaller chemical speed, 
which may be surprising since the chemical speed, con- 
trolled by a diffusive flux, is expected to be larger when 
the gradients are sharper. Actually, the chemical is not 
larger than in the sharp gradient case, precisely because 
of the difference in the initial gradients magnitude, but 
it increases much faster. 



IV. CONCLUDING REMARKS 

We have studied an infinitely fast bimolecular chemi- 
cal reaction in a two-dimensional Navier-Stokes flow at 
moderate Reynolds number with chaotic advection. The 
computation of the probability distribution function of 
the Lyapunov exponents suggests that large deviation 
theories may be relevant to describe its behavior after 
a few turnover times. We defined G e (X, t) such that the 
FTLE pdf scales like e' tG ^ x ^ and minG e (A,t) = 0. The 

function G e satisfactorily converges to a Cramer function 
G in a couple of turnover times, at least for exponents 
larger than their mean value. 

We have investigated the early regime (~ 5 turnover 
times of the flow) of the reaction, corresponding to the 
time window where the contact line is a clearly defined 
material line that does not depend on diffusion. We pos- 
tulate that this time window is limited by the mix-down 
time scale from the large scales to the diffusive cutoff 
and scales like the log of the Peclet number. We have re- 
lated, both theoretically and numerically, the Lagrangian 
straining properties of the flow, as captured by the joint 
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Time in units of turnover time T 

FIG. 10. (Color online) Ensemble average of the chemical 
speed, in the smooth gradient case, divided by k. The sym- 
bols correspond to numerical results from the 34 members en- 
semble, for different Prandtl numbers Pr — 2, 4, 8, 128. The 
solid line (calculation from the FTLE) correspond to (|42[) . 
The exponential increase at a rate 0.09 corresponds to the 
expected asymptotic regime of (|42[1 . as expressed in (|43[) and 
has been shifted vertically for clarity. We note the log scale 
in the y-axis. 



pdf of the Lyapunov exponents A and two equivalent 
times r and r (f30j) . to the following quantities: 

• The ensemble average contact line length be- 
tween the reactants (L). After a brief tran- 
sient corresponding to the equilibration of the con- 
tact line with the flow, i.e. to the alignment of 
the contact line elements with the direction cor- 
responding to the maximum growth, independent 
of its initial orientation, the contact line lengthens 

max[A-G e (A,t)]t . , 

like e x which converges m time to e 1 , 

where Ai is the Legendre transform of G evaluated 
in one and is determined by rare large events in the 
FTLE distribution. 

• The ensemble mean of the gradients along 
the contact line (|V(/>£|). It scales like k" and is 
determined by the pdf of (A, r, f ) through ([57)1 . The 
influence of r diminishes with time as the contact 
line is equilibrating with the flow. The dependence 
between A and r is crucial to accurately predict 
(|V(/)£|). Our main assumption was the stationar- 
ity of the Lyapunov vectors, justified by their fast 
exponential convergence in time. It would be inter- 
esting to extend this work without this assumption 
to precise the conditions of its applicability. 



• The ensemble mean chemical speed. The 

chemical speed is defined as the modulus of the time 
derivative of the sum of the two reactants' mean do- 
main concentrations. It scales like in the limit 
of infinite initial gradients. This scaling is consis- 
tent with Q in the special case of a contact line 
of dimension one separating two on/off fields. The 
ensemble average chemical speed is proportional to 
the product of (L) and (\V(f>c\}- Hence, an initial 
decrease of the chemical speed is related to the de- 
crease of the gradients, while a later regime is dom- 
inated by the lengthening of the contact line and 

. . max[A-G c (A,t)]t 

is consequently equivalent to e A . Both 

the contact line length and the chemical speed are 
determined by very rare events in the tail of the 
FTLE distribution. This points out the importance 
of considering the distribution of the FTLE, which 
is not al way s taken into account in the literature 

The case of smooth gradients exhibits some significant 
differences. The gradients increase instead of decreasing 
and are initially not affected by diffusion. The two main 
consequences are that the chemistry scales like n and 
increases exponentially in time at a rate determined by 
even rarer events in the tail of the FTLE distribution 
63). 

The theory developed in this paper should allow to 
predict the evolution of the pdfs of the gradients along 
the contact line and of the passive tracer which would 
be a very robust way to test it. This will be the subject of 
a future paper. Another paper in preparation which takes 
into account the fractal structure of the contact line in 
order to investigate the intermediate regime, where the 
chemical production reaches a maximum, and the long 
term decay of the reactants. 

Some interesting open questions about the Lagrangian 
properties of a two-dimensional Navier-Stokcs flow have 
arisen from this study. What determines the initial time 
evolution of the FTLE pdf? What determines the shape 
of the Cramer function Gl Is it possible to predict the 
asymptotic form of the pdf of -? More importantly, the 
dependence between r and A, associated with the persis- 
tence of significant probability of small FTLE, seems to 
be a major difference with simple prescribed flows used in 
the literature to study chaotic advection, which may not 
be, as a consequence, representative of dynamically con- 
sistent flows. Studying the joint pdf of (A, r) may happen 
to be useful to better understand the mixing of both pas- 
sive and active tracers in two-dimensional Navier-Stokcs 
flows with chaotic advection, particularly when using La- 
grangian straining theory approaches. 
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